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The  impact  of  a symmetrical  array  geometry,  the  use  of  a quantized  stored,  cosine  function,  the 
exploitation  of  digital  Fourier  transform  algorithms,  and  the  application  of  trigonometric  interpolation  in 
the  computation  of  array  patterns  is  discussed.  Careful  selection  of  parameters  permits  sampling  the  array 
pattern  only  6%  above  the  theoretical  Nyquist  limit.  Reconstruction  of  array  patterns  showing  —20,  —30, 
and  — 40-dB  relative  interpolation  errors  are  presented.  A saving  of  8000:1  in  computation  time  over  direct 
“brute-force”  array-pattern  computation  is  illustrated  for  a hypothetical  array. 

PACS  numbers:  43.60.Qv,  43.60.Cg,  43.30.Vh 


INTRODUCTION 

In  today’s  world  of  readily  available  high-speed 
digital- computer  facilities,  the  computation  of  the  di- 
rectional response  for  arrays  composed  of  many 
arbitrarily  distributed  elements  is  commonplace.  How- 
ever, when  the  number  of  elements  is  very  large  and 
the  aperture  extends  over  many  wavelengths  in  two  or 
three  dimensions,  the  cost  of  such  computation  can  be 
significant.  When  the  degree  of  array  complexity  is 
reached  where  computing  cost  becomes  a concern,  a 
careful  selection  of  algorithms  and  sampling  intervals 
for  numerical  processing  can  lead  to  substantial  sav- 
ings in  computer  time. 

The  rationale  for  the  selection  of  sampling  intervals 
begins  with  the  basic  transform  relationship  embodied 
in  the  farfield  directional-response,  function  of  a co- 
herently excited  aperture; 

4>^)=  [ >l(p)exp[f27ric -pidp  . ^ (1) 

When  the  aperture  is  made  up  of  N discrete  phased 
array,  elements,  A(^  consists  of  a set  of  delta  func- 
tions with  complex  amplitude  coefficients.  For  this 
case,  which  is  the  one  of  interest  here, 

N 

«5(k)=  A.extpiiltik-pj]  . (2) 

i=i  • - . 

Furthermore,  if  the  array  is  phased  for  plane  waves, 
the  complex  phase  of  the  A^  coefficients  can  be  in- 
corporated into  k so  that,  in  Cartesian  coordinates, 
the  array  factor^  is  defined  as 

N 

4>^)=  S rtijexp[t2Tr(k^Xj+k,,yj+k^Zj)]  , (3) 

i=i 

where  {a  - ao  )/X;  k^=  {13-  /3o)A;  k^=  (y_  \ is 

the  wavelength  of  the  plane  wave;  a,  (3,  y are  the  direc- 
tion cosines  of  the  normal  to  the  wavefront;  and  a^, 

Po,  To  are  the  direction  cosines  of  the  steering  vector. 
Note  that,  in  accord  with  conventional  antenna  terminol- 
ogy, a distinction  is  made  between  the  generally  com- 
plex array  factor  tpik)  and  its  power,  the  array  pat- 
tern l;/)(k)l^. 

The  array  factor  provides  a convenient  way  of  de- 
scribing the  performance  of  a distributed  set  of  trans- 
mitting or  receiving  elements.'  However,  it,  of  course, 
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does  not  include  the  directional  characteristics  of  the 
elements  themselves.  On  the  one  hand,  if  all  of  the, 
elements  are  identical  and  oriented  with  their  response 
axes  parallel,  the  overall  response  factor  at  a particular 
frequency  is  a simple  product  of  the  array  factor  and 
the  directional  factor  of  the  individual  elements.  On 
the  other  hand,  if  they  are  not  identical  or  not  oriented 
parallel  to  each  other,  the  element  weights  are  not  in- 
dependent of  angle.  In  this  case  the  computation  be- 
comes much  more  laborious  in  that  the  response  is  no 
longer  a function  of  the  universal  parameters  k^,  ky, 
kg  but  must  be  considered  a function  of  frequency,  re- 
sponse direction,  and  steering  direction. 

The  array  factor  ^(k)  is  a continuous  function  of  k 
but,  in  practice,  it  is  always  approximated  by  a dis- 
crete set  of  M samples  over  the  “visible”  range  of  1? 
that  is  of  interest.  (The  visible  range  is  that  range  of  . 
values  corresponding  to  sums  of  the  direction  cosines 
for  the  response  and  steering  vectors  lying  between 
4-  2 and  - 2 at  the  frequency  of  interest.  Thus,  the 
straightforward  computation  of  Eq.  (3)  requires  A/xN 
operations  of  the  form 

ARC  = KX(J)  * X(I)  K Y(J)  ♦ Y(1  -F  KZ(J)  ♦ Z (1) , 

CSUM(J)  = CSUM(J)-f  A(l)*  COS(ARG)  , (4) 

SSUM(J)  = SSUM(J)-fA(1)*S1N(ARG)  . 

I.  EFFECT  OF  SYMMETRY 

For  the  sake  of  completeness,  what  is  perhaps  the 
most  common  technique  for  reducing  the  computation 
effort  will  be  mentioned— that  of  imposing  a condition 
of  symmetry  on  the  array.  Consider  an  array  distri- 
bution which  is  symmetrical  about  all  three  axes.  In 
this  case  there  will  be  eight  elements  described  by  one 
set  of  coordinate  values  (±  \x\,  ± lyl,  ± Iz  l).  Suitably 
grouping  the  terms,  the  eight-element  sum  becomes 

^exp(±ffe,j|x|  ±ffey(y|  ±ikg\z\) 

e 

= 8cos(fe,j|x|)cos(fey|y|  )cos(fe,,|2|)  . (5) 

The  eight  complex  terms  of  the  sum  thus  combine  into 
a single  symmetrical  term  such  that  the  imaginary  part 
of  the  sum  equals  zero  and  need  not  be  computed. 
Furthermore,  because  is  now  symmetrical,  i.e.. 
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kyf  I J I I ^ J 

just  its  real  part,  and  that  for  only  positive  values  of 
k^,  ky,  need  be  computed.  If  the  array  is  not  sym- 
metric, only  the  redundancy  :/)  (k  ) = :/)*(- k)  holds  so  that 
the  complex  ip^)  must  be  computed  over  one-half  of  the 
appropriate  k space.  Thus,  for  the  symmetrical  array, 
the  2xM/2xn  operations  for  a three-dimensional  array 
and  a three-coordinate  beam  set  of  the  form  of  Eq.  (4) 
have  been  replaced  by  (M/8)x  (N/8)  operations  of  the 
form 

SUM(J)  = SUM(J)  + A(I)  ♦ COS[KX(J)  * X(I)] 

♦ C0S[KY(J)*  Y(I)]+COS[KZ(J)*  Z(D]  . (6) 

Recognizing  that  the  COS-function  routine  dominates 
the  computational  time,  the  overall  reduction  in  time 
for  a three-dimensional  array  is  about  ^ . For  two- 
and  one-dimensional  arrays  the  saving  is  not  quite  as 
great,  about  i and  i,  respectively.  Without  question, 
this  order-of-magnitude  savings  in  time  is  a strong 
argument  for  symmetry  in  array  design  as  far  as  com- 
putation of  array  patterns  is  concerned. 

II.  USE  OF  A STORED  COSINE  TABLE 

Another  very  effective  way  of  reducing  the  computa- 
tion time  is  to  replace  the  COS-function  routine  with 
a look-up  table.  This  inevitably  will  introduce  an 
error  in  the  result  because  of  the  quantization  intro- 
duced. The  statistical  effect  of  the  quantization  error 
can  be  estimated  by  using  the  expression  for  the  effect 
of  random  phase  errors  in  an  array,  ^ 

Jil 

where  i/)o(fe)  is  the  array  pattern  without  errors,  is 
the  variance  of  the  phase  ^errors,  <r\  is  the  variance  of 
the  amplitude  errors,  and  the  are  the  amplitude  co- 
efficients of  the  array  elements.  An  insight  into  the 
magnitude  of  the  effect  of  phase  errors  can  be  gained 
by  considering  an  array  with  elements  having  unit 
weights.  The  difference  between  the  mean  pattern  and 
the  pattern  unperturbed  by  errors  will  be  given  by 

f{k)~  4>l (k)  = crl [N -4>l(k)]  + a\N.  (8) 

On  axis  (k=0),  ipo{k)  = N^  so  that  represents  the  frac- 
tional loss  in  power  in  the  main  lobe. 

If  the  elements  of  the  array  were  randomly  dis- 
tributed, the  mean  side-lobe  level  [ipoik)]  would  be 
equal  to  N.  The  average  effect  of  phase  errors  which 
further  randomize  the  phase  of  the  element  summations 
is  essentially  zero  for  this  case.  In  regions  of  side 
lobes  where  ipl(k)  is  larger  than  N,  the  effect  of  phase 
errors  will  be  to  reduce  the  side-lobe  level;  where  ipoik) 
is  less  than  N,  phase  errors  will  increase  the  side- 
lobe  level.  A useful  measure  of  the  quantization  error 
is  obtained  by  placing  ipQ{k)  = 0,  then  the  residual  mean 
error  pattern  is  a^N  and  cr^  represents  the  fractional 
quantization  noise  power  in  the  side-lobe  region  as  ref- 
erenced to  the  mean  side-lobe  power  level.  If  A<^  is 
the  step  size  in  the  argument  quantization,  the  error  can 
be  treated  as  being  uniformly  distributed  over  the  inter- 
val (-2A0-+|a0).  For  this  distribution,  cr|  = -R  A0^. 
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Thus,  for  a quantization  noise  that  is  20  dB  below  the 
mean  side-lobe  level,  = 0.  01,  A0  0.  35  rad.  Fora 
level  30  dB  below  the  mean  side-lobe  level,  A0  ^0. 1 
rad.  It  is  apparent  that,  for  most  purposes,  a rela- 
tively crude  quantization  of  the  look-up  table  may  be 
used— intervals  between  0.  1 and  6.  3 rad.  Depending 
upon  the  architecture  of  the  specific  computer  in- 
volved, it  may  be  more  advantageous 'to  use  either  a 
short  look-up  table  that  covers  the  range  of  27t  in  2" 
intervals  and  use  the  lowest  n bits  of  the  argument  for 
the  address,  or  else  use  an  extended  table  covering  the 
entire  range  of  arguments,  i.e.,  2ttD/X,  where  D is 
the  maximum  dimension  of  the  array  aperture.  In 
either  case,  the  savings  in  time  will  result  from  re- 
placing the  operations:  three  cosine  routines  plus  six 
multiplications  plus  one  addition,  by  the  operations: 
six  multiplications  plus  one  addition.  If  the  cosine 
routine  takes  typically  ten  times  as  long  as  the  multi- 
plication routine,  a reduction  in  processing  time  of 
about  I will  be  achieved  by  the  use  of  a look-up  table. 


III.  UNIFORM  SAMPLING  OF  THE  ARRAY  PATTERN 

From  Eq.  (1)  it  is  clear  that  ipO^)  and  A(p)  are  three- 
dimensional  Fourier-transform  pairs.  In  order  to 
treat  the  sampling  of  the  variables  separately,  the 
volume  integration  of  Eq.  (1)  can  be  rewritten 


a^exp{i2trk^x)dx  , 


where 

a^=a^(x,ky, 


k^)=  { f A{x,  y,  z) 


(9) 


Xe-xp[i2-n{kyy  + k^z)]dy  dz 

and  fly  and  can  be  similarly  defined.  The  resulting  ■ 
one-dimensional  transform  pairs  are  of  the  form 


tj){k)=  I a,,exp{i2Ttk^x)  dx  , 

— ce 

r ip(k)exp{- i2irk^x)dx  . 


(10) 


It  is  apparent  that  the  physical-aperture  function 
is  the  Fourier  spectrum  of  the  array  factor  i>{k).  If 
ip{k)  is  to  be  sampled  at  uniform  intervals  of  Ak,  the 
resulting  sampled  function  can  be  expressed  as 

n=* 

i>^(k)  = ip(k)  Yj  b{k-nAk)  . (11) 


Here  the  sampled  function  is  obtained  by  multiplying 
the  continuous  function  by  a Dirac  comb.  * The  prop- 
erties of  the  Dirac  comb  are  such  that  the  spectrum  of 
the  sampled  function  will  be  given  by 


asx^)=  Y axi^-  . (12) 

The  rbsult  is  the  familiar  aliasing  or  replicating  of  the 
original  spectrum  at  frequency  intervals  in  x of  {Ak)~^, 
Because  the  physical  aperture  is  clearly  restricted  to 
lie  within  its  maximum  dimensions,  the  maximum  ex- 
tent of  the  spectrum  of  ip(k)  is  (~^D^x^\d).  Thus, 
if  Ak^\/D  (Nyquist  period)  the  aliased  spectra  will  not 
overlap. 
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It  is  customary  to  reconstruct  tp(k)  from  i>^{k)  by 
selecting  the  baseband  (m  = 0)  portion  of  the  aliased 
spectrum  by  multiplying  by  a filter  window  GOf).  Then 
the  reconstructed  i>(k), 

4>R{k)=  f G(x)a^^{x)e7i^(-i2i!kx)dx  (13) 

J 

or  in  its  equivalent  convolution  expressions, 


-C)d^ 


-<sc 


4>R(k)=  H{k-n^k)  , 


(14) 


. where 


H(k)  = 

1 Gi^)e'xp(—i2-nkx)dx, 

/ - co  ^ 

(15) 

G(r)  = 

f H(k)e-xp(i2-nkx)dk  . 

/ «eO 

(16) 

As  shown  by  Eq.  (14),  the  reconstruction  operation 
is  one  of  interpolation;  that  is,  the  value  of  ip  for  a 
value  of  k lying  between  n^k  and  (w+  .1)  Afe  is  obtained  as 
a weighted  sum  of  the  samples  of  ip  around  that  value. 

The  most  desirable  shape  for  G{x)  would,  of  course, 
be  a rectangular  window 

g(at)=i  (in^x^in), 

= 0 (iD>x>iB). 

Unfortunately,  as  is  well  known,  the  /I(k)  function  for 
such  a window  (sm-nkD)/-nkD  extends  over  a range  of 
i"®.  In  order  to  carry  out  the  reconstruction  in  apracti- 
cal  way,  it  is  necessary  to  restrict  the  interpolation 
sum  of  Eq.  (14)  to  a limited  set  of  samples  around  the 
interpolation  point.  This  will  give  rise  to  an  imperfect 
reconstruction  because  the  aliased  spectrum  cannot  be 
totally  rejected.  As  a measure  of  the  imperfection,, 
the  rms  deviation  of  the  filtered,  aliased  spectrum  from 
the  original  spectrum  will  be  computed  as  follows:  As- 
suming a uniformly  weighted  aperture,  ajx)=  1 and  a 
G(x)  normalized  to  a mean  value  of  1.  0 at  a:  = 0, 

ERR=  i^^^'\G{x)-\fdx+  G^(x)dx^  j . 

(17) 

This  error  serves  as  a criterion  for  the  selection  of 
both  a suitable  sampling  period  A)fe  and  an  appropriate 
interpolation  algorithm. 


IV.  SOME  INTERPOLATION  ALGORITHMS 

Several  methods  of  interpolation  are  of  interest  for 
the  reconstruction  process.  These,  along  with  their 
associated  G{x)  spectral  windows  and  their  interpola- 
tion errors  are  given  below. 

A.  Boxcar  smoothing 

Whenever  the  continuous  function  :/)(fe)  is  represented 
by  a set  of  discrete  values,  this  form  of  interpolation 
is  imposed.  That  is,  the  function  is  considered  to  have 
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the  value  in  the  interval  ikj—  k^k<k<kj+jAk).  The 
algorithm  is  defined  as  , 

/f^(/fe)=l,  0«  l^l  ^jAk  , 

H^{k)=0,  \Ak<\k\  . (18) 

Its  spectral  window  is 

G = sin{-nx Ak) / -nx Ak  . 

The  magnitude  of  this  window  function,  shown  in 
Fig.  1,  features  an  infinite  set  of  nulls,  each  centered 
on  a multiple  of  the  sampling  frequency.  Thus  for  a 
sufficiently  high  sampling  rate,  the  relative  width  of 
the  aliased  spectra  of  ip,  which  are  centered  around 
multiples  of  the  sampling  frequency,  can  be  made 
small  enough  to  achieve  any  desired  degree  of  rejec- 
tion. This  is  the  ultimate  interpolation  imposed  on  the 
sampled  function  after  a sufficiently  high  effective 
sampling  frequency  has  been  achieved  by  more  efficient 
interpolation  algorithms.  The  errors  inherent  in  the 
gradual  roll  off  within  the  desired  pass  band  are  cor- 
respondingly reduced  by  an  increase  in  effective  sam- 
pling frequency. 


0 0.5  1.0  1.5  2.0  2.5  3.0  3.5  4.0 


Frequency  X.  1/A  K 

FIG.  1.  (a)  Frequency  window  G(x)  for  boxcar  interpolation 
and  (b)  frequency  window  Gix)  for  linear  interpolation. 
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B.  Linear  interpolation 

In  this  method  the  values  of  the  function  at  each  end 
of  a sampling  interval  are  combined  with  proportional 
weights  to  obtain  intermediate  values.  The  common 
application  of  this  method  is  the  practice  of  connecting 
points  on  a graph  by  straight-line  segments.  The  win- 
dow functions  are 

Hg{k)=l-\k/i^k\  (0<  l/fe|  « A/fe)  , 

= 0 (A/fe<  \k\  , (19) 

Gg(x)  = [sin(Trx^k)/Trx^kf  . 

This  window,  also  plotted  in  Fig.  1,  exhibits  deeper 
nulls,  at  the  same  multiples  of  the  sampling  frequency, 
than  the  previous  one.  It  suffers,  however,  from  a 
greater  distortion  within  the  pass  band  inasmuch  as  the 
roll  off  in  that  region  is  also  greater  than  that  for  the 
boxcar  smoothing. 

C.  Trigonometric  interpolation 

The  trigonometric- interpolation  method  is  ideally 
suited  to  uniformly  spaced  data  samples  and  is  closely 
related  to  the  Fourier-series  expansion.  The  method 
is  described  by  Lanczos®  in  some  detail.  His  formula- 
tion can  be  used  to  derive  an  expression  for  the  weight- 
ing coefficients  corresponding  to  the  window  function 

In  brief:  Given  a set  of  n+  1 uniformly  spaced  sam- 
ples fj  of  an  arbitrary  function  f(pc),  first  scale  x such 
that  the  interval  O^x^n  corresponds  to  the  sampling 
interval  0 « 7 « M.  Now,  generate  a new  function 

h(^)  = f(x)-Ci.-x/n)fa-ipc/n)f„.  (O^x^n),  (20) 

h(-x)  = -h{x)  i-n^x^O). 

hM  is  now  an  odd  function  of  x.  It  is  identically  zero 
at  —n,  0,  n and  the  first  derivative  is  continuous  across 
these  nodes.  Thus,  it  can  be  readily  expanded  in  the 
pure  sine  series 

*(^)=  6^sin(^']  , (21) 

^=1  \ n / 

where 


TABLE  I.  Trigonometric  weights  for  midinterval  interpola- 
tion.^ 


m 

2 

k + \) 

4 

8 

16 

1 

0.500 

0.5773 

0.6259 

0.  6343 

2 

-0.0773 

-0.1791 

-0.2052 

3 

0.0688 

0. 1155 

4 

-0.0156 

-0.0740 

5 

0. 0484 

6 

-0.0297 

7 

0.0142 

8 

-0.0035 

*In  this  table  («  + l)  is  the  number  of  samples  used  for  interpola- 
tion and  m is  the  index  of  symmetrical  displacement  from  the 
central  interval. 
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Combining  (20) -.(22) 


Xsin(f^)|  + /o(l-^)+^/„.  (23) 

We  are  interested  in  a set  of  interpolation  weights  such 
that 

n 

fM=Y^(P,{x)f,  (24) 

g=0 

collecting  coefficients  in  Eq.  (23)  the  set  of  w + 1 weights 
can  be  written  as 


,(«)=- X]  sin  sin  (l«?«w-l), 

n i„i  \n  J \n  / 

n \n  ) \n  ) 


It  is  interesting  to  note  that,  if  only  two  samples  are 
used,  n-  \ and  the  summations  have  zero  terms  so  that 
the  interpolation  reverts  to  /(»:)=  x)fQ  + xfi,  which  is 

just  the  linear  interpolation  mentioned  above.  Fur- 
ther, by  virtue  of  the  orthogonality  of  the  functions,  if 
X is  an  integer  corresponding  to  one  of  the  sample 
points,  the  weight  for  that  sample  point  equals  one  and 
all  the  other  weights  are  equal  to  zero,  confirming  that, 
indeed,  the  expansion  reproduces  the  sample  points 
exactly. 


For  any  value  of  x within  a sampling  interval,  a set 
of  (p’s  may  be  calculated.  In  our  application  of  the 
trigonometric  interpolation  the  weighting  coefficients 
for  the  middle  of  the  sampling  interval  will  be  used. 
These  weights  are  given  in  Table  I for  interpolation 
windows  utilizing  2,  4,  8,  arid  16  samples. 


The  addition  of  a set  of  midinterval  interpolated  points 
doubles  the  effective  sampling  frequency.  This  in  ef- 
fect aliases  the  continuous  G(.x)  window  associated  with 
the  continuously  distributed  weights  of  Eq.  (25)  around 
multiples  of  2/A,%,  twice  the  original  sampling  frequency. 
The  aliased  windows  can  be  computed  empirically  from 
an  FFT  transform  of  the  delta-function  sequence  made 
up  of  the  set  of  midinterval  weights  superimposed  on  the 
integer  interval  weights,  e.g. , for  (w  + l)  = 4 the  set  of 
uniformly  spaced  delta  functions  to  be  transformed  have 
amplitudes  of  -0.0773,  0,  0.5773,  1.0,  0.5773,  0, 
-0.0773.  These  aliased  windows  are  shown  in  Fig.  2 
for  the  set  of  interpolation  weights  of  Table  I.  In  each 
case  the  aliased  pass  bands  centered  at  2n/^k  are  ap- 
parent. As  one  would  expect,  increasing  the  number 
of  samples  improves  (flattens)  the  acceptance  band  and 
also  improves  (steepens)  the  rejection  band. 


For  reasons  stated  later,  the  complete  interpolation 
will  be  performed  as  a sequence  of  midinterval  inter- 
polations. The  first  step  of  the  sequence  effectively 
doubles  the  sampling  frequency,  rejecting  the  odd 
aliased  aperture  bands  but  retaining  the  even  ones  in 
accord  with  the  window  functions  of  Fig.  2.  The  set  of 
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Frequency  X,  1/AK 


Frequency  X,  1/AK 


FIG.  2.  (a)  Aliased  frequency  window  for  2-weight  trigono- 
metric interpolation,  (b)  aliased  frequency  window  for  4-weight 
trigonometric  interpolation,  (c)  aliased  frequency  window  for 
8-weight  trigonometric  interpolation,  and  (d)  aliased  frequency 
window  for  16-weight  trigonometric  interpolation. 


y sampling  frequency  ratio 

FIG.  3.  Relative  rms  interpolation  errors  versus  the  ratio  of 
the  sampling  frequency  to  Nyquist  (D).  Errors  observed  for 
the  reconstructions  of  Figs.  4,  5,  and  6 are  plotted  as  circles 
for  comparison. 

samples  now  represent  aperture  bands  aliased  at  2w/ ^k, 
but  the  second  interpolation  step  will  be  operating  at  an 
effective  sampling  frequency  of  4/ tik  for  m odd  and  re- 
ject the  bands  for  m even.  Thus,  the  stop  bands  of 
successive  midinterval  interpolation  windows  are  pre- 
cisely aligned  with  the  pass  bands  of  the  previous  step. 
Also,  of  course,  the  relative  width  of  the  aperture  band 
with  respect  to  the  sampling  frequency  is  reduced  by  a 
factor  of  2 with  each  iteration  step  so  that  lower -order 
windows  may  be  used.  The  final  iteration  step  in  every 
case  will  use  one  of  the  windows  of  Fig.  1,  i.e.,  either 
the  boxcar  window  which  corresponds  to  the  representa- 
tion of  the  array  factor  or  of  the  array  pattern  by  a set 
of  numbers,  or  the  linear  window  which  corresponds  to 
the  representation  of  the  pattern  by  a continuous  plot 
with  straight-line  segments  connecting  the  adjacent 
interpolated  points.  Each  of  these  final  windows  has  a 
rejection  null  centered  on  all  of  the  residual  aliased 
aperture  bands  resulting  from  the  prior  succession  of 
interpolation  steps. 

V.  INTERPOLATION  ERRORS 

The  precision  of  the  interpolation  restoration  can  be 
estimated  by  numerically  integrating  Eq.  (17)  over  the 
response  curves  of  Figs.  1 and  2.  These  errors  are 
presented  in  Fig.  3 as  a function  of  y,  the  ratio  of  the 
effective  sampling  frequency  to  Nyquist;  i.  e. , units  of 
D.  Here  one  sees  the  dramatic  results  of  the  trigono- 
metric interpolation. 

As  an  aid  to  formulating  a complete  interpolation 
algorithm,  the  intersections  of  the  error  curves  of  Fig. 
3 with  various  error  levels  are  tabulated  in  Table  H. 

As  the  number  of  interpolation  points  is  increased  there 
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TABLE  n.  Minimum  sampling  frequency  ratio  y for  interpola- 
tion errors. 


Number  of  points  used  for  interpolation 

Root-mean-square 


error  (dB) 

1 

2 

4 

8 

16 

-20 

4.62 

2.08 

1.35 

1. 11 

1.01 

-30 

13.80 

3.63 

1.78 

1.23 

1.06 

-40 

39.80 

6.31 

■ 2. 19 

1.74 

1.26 

must  be  a corresxxDnding  addition  to  the  number  of  data 
points  computed.  Thus,  if  the  number  of  samples  of  k 
which  are  to  be  computed  is  small  enough,  the  advan- 
tage gained  in  the  reduction  of  sampling  frequency  by 
using  more  points  in  interpolation  is  offset  by  the  num- 
ber of  additional  points  required  for  the  interpolation 
filter.  The  minimum  ranges  of  k in  Nyquist  intervals 
(1/Zl)  at  which  a given  interpolation  filter  will  save 
computation- time  are  given  in  Table  III  for  the  condi- 
tions of  Table  II.  These  are  limiting  values  obtained  by 
assuniing  that  the  number  of  elements  in  the  array  is 
very  large  compared  to  the  number  of  points  used  for 
interpolation.  Except  for  very  crude  and  broad  pat- 
terns (-20  dB  error,  -5  <kD<5)  interpolation  with  8 or 
even  16  weights  is  worthwhile.  A 16-weight  window  will 
allow  sampling  at  1. 06Nyquistwith  anerrorof  - 30  dB. 

VI.  INTERPOLATION  SEQUENCES 

As  an  .example  of  the  selection  of  a complete  inter- 
polation algorithm,  consider  a requirement  for  an  er- 
ror of  -40  dB  and  a range  for  kD  of  20.  From  Table 
ni  we  select  a 16 -point  filter  for  the  first  step,  from 
Table  n the  frequency  ratio  y=  1. 26.  The  first  step  inr 
creases  y to  2.  52.  Again  from  Table  II  the  next  step 
can  use  a 4-point  filter  to  increase  y to  5.04.  How- 
ever, this  frequency  ratio  is  below  that  required  for  the 
2-point  interpolation  so  the  4-point  filter  has  to  be  re- 
peated to  give  a y of  10.08.  At  this  point  the  array -fac- 
tor data  can  be  plotted  since  the  2-point  filter  is  in- 
herent in  the  plotting  process.  , If,  on'the  other  hand, 
one  wishes  to  represent  the  array  factor  as  a numerical 
data  array,  y has  to  be  greater  than  39.  8 to  retain  the 
-40-dB  error  band.  Thus  the  sequences  are  as  follows: 

-40  dB  (plotting) : 16,  4,  4, 

- 40  dB  (data  array) : 16,  4,  4,  2,  2 . 

This  example  brings  out  two  reasons  why  the  method 
of  sequential  midinterval  interpolation  steps  was  chosen 
rather  than  some  other  combination  of  fractional  in- 
terval interpolation.  First,  a single,  midinterval  in- 
terpolation step  brings  the  process  well  within  range 
for  a lower -order  simpler  window.  Second,  the  mid- 
interval 2 sample  window  for  a final  iteration  step  is 
much  simpler  to  implement  than  any  other  fractional 
interval  in  that,  with  integer  arithmetic,  it  consists  of 
a simple  binary  sum  and  right  shift.  This  can  be  quite 
important  because  this  last  step  consumes  a signifi- 
cant portion  of  the  interpolation  computing  time  as  will 
be  evident  in  the  example  later. 

Before  illustrating  the  interpolation  process,  a point 


should  be  raised  regarding  the  spectrum  of  the  array 
pattern.  Since  the  array  pattern  is  derived  from  the 
square  of  the  array  factor  it  has  a spectral  width  twice  ' 
as  great  as  the  aperture  D.  Thus,  the  array  pattern 
requires  a sampling  frequency  twice  that  for  the  array 
factor  in  order  to  avoid  aliasing.  In  order  to  accommo- 
date the  broader  spectrum  an  additional  step  of  mid- 
interval interpolation  should  be  imposed  before  con- 
verting the  array  factor  to  an  array  pattern. 

VII.  PLOTTING  EXAMPLE 

An  interpolation  sequence  of  8,  4,  2 has  been  chosen 
for  the  purpose  of  illustration.  Figures  4,  5,  and  6 are 
[sin(7T^)/77^]*  array  patterns  interpolated  by  the  use  of 
the  above  sequence  on  the  array  factor,  siniw k)/ it k.  The 
squares  (plotted  with  a +10-dB  offset)  represent  the 
data  points  used  for  interpolation  and  the  line  connect- 
ing these  points  is  the  uninterpolated  plot.  Sampling 
frequency  ratios,  referred  to  D,  the  Nyquist  frequency, 
of  1.1  (Fig.  4),  1.25  (Fig  5),  and  1.75  (Fig.  6)  were 
used.  The  rms  deviations  of  the  interpolated  data  points 
from  the  true  values  are  superimposed  on  the  8 -point 
curve  of  Fig.  3.  These  show  good  agreement  between 
the  predicted  errors  and  those  actually  realized.  The 
three  cases  correspond  to,  roughly,  -20-,  -30-, 
and  —40-dB  errors. 

Although  the  central  lobe  of  the  interpolated  array  pat- 
tern of  Fig.  4 (1.  IX  Nyquist)  matches  the  [sin{kx) / kxf 
pattern  quite  closely,  the  deviations  of  the  side-lobe 
peaks  from  the  theoretical  values  are  quite  large.  In 
the  region  of  the  dip  in  the  side -lobe  pattern  an  error 
of  6—7  dB. occurs.  This  is  not  inconsistent  with  the 
-20-dB  error  level  because  the  side-lobe  level  itself 
is  down  to  — 30  dB  and  only  small  absolute  errors  are 
required  to  change  the  relative  levels  of  these  low  side 
lobes  significantly. 

At  the  slightly  higher  sampling  frequency  of  Fig.  5 
(1.  25  Nyquist)  the  side-lobe  behavior  is  much  better, 
with  errors  on  the  order  of  1 dB.  The  side -lobe  levels 
of  Fig.  6 (1.75  Nyquist)  are  correct  to  within  about  0.2 
dB. 

VIII.  UNIFORM  SAMPLING  OF  THE  ARRAY 
APERTURE 

At  this  point  it  is  evident  that  both  the  aperture  and 
array-factor  fields  are  represented  by  sets  of  discrete 
data  samples.  The  array-factor  set  is  generally  uni- 
formly sampled  but  the  aperture  set  is  not.  If  both 
sets  were  to  be  composed  of  the  same  number  of  uni- 
formly spaced  samples,  the  transform  of  Eq.  (1)  would 

TABLE  in.  Minimum  fe/Afe  range  (Nyquist  intervals)  at  which 
a given  interpolation  filter  saves  computation  time. 

Number  of  points  used  for  interpolation 

Root-mean-square 


Error  (dB) 

1 2 

4 

8 

16 

-20 

0.8 

3 

15 

80 

-30 

0.2 

i.  1 

7.3 

47 

-40 

0.06 

0.5 

9 

17 
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FIG,  4.  Array  pattern  generated  at  a sam- 
pling frequency  of  1. 1 x Nyquist.  Squares 
are  the  computed  pattern  points  and  the  as- 
sociated uninterpolated  plots  are  offset  +10  dB. 
The  solid  curve  is  derived  by  an  8,  4,  2 
trigonometric-interpolation  sequence  frorri 
the  same  set  of  computed  data  points 
(squares).  The  relative  rms  error  is  —20.6 
dB. 


take  on  the  form  of  a discrete  Fourier  transform  (DFT) 
and  the  powerful  fast-Fourier -transform  (FFT) 
algorithms  could  be  exploited  in  the  computation.  As 
the  aperture  size  and  the  number  of  elements  are  in- 
creased the  computational  efficiency  of  the  FFT  al- 
gorithms becomes  more  important. 

One  approach  to  obtaining  a uniformly  sampled  array 
aperture  is  to  place  the  elements,  by  design,  on  the 
interstices  of  a coordinate  lattice  corresponding  to  an 
appropriate  uniform  sample  set.  This  approach,  al- 
though simple  and  effective,  is  unduly  restrictive  in  the 
selection  of  element  positions.  As  an  alternative, 
for  any  arbitrary  array,  each  delta  function  of  the  ar- 


ray-element distribution  must  be  expanded  into  a set 
of  uniformly  spaced  spatial  samples  in  the  vicinity  of 
the  element.  These  individual  sets  are  then  summed  to 
obtain  the  total  array-aperture  distribution.  This  ex- 
pansion is,  in  essence,  the  inverse  of  the  interpolation 
of  the  previous  section  where  the  value  at  an  arbitrary 
point  was  obtained  from  the  values  of  neighboring  ele- 
ments of  a uniformly  sampled  data  set. 

Because  the  value  of  the  delta  function  is  zero  at  both 
ends  of  the  expansion  interval,  a trigonometric  sine- 
series  expression  is  appropriate.  The  approach  is  to 
transform  the  delta  function  into  a finite  set  of  sine- 
series  coefficients  and  then,  through  an  inverse  dis- 
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Kappa,  A K 


FIG.  6.  The  same  presentaKon  as  Fig. 
4 at  a sampling  frequency  of  1. 75  x Ny- 
quist.  The  relative  rms  error  is  -41.5 
dB. 


Crete  Fourier  Transform  obtain  a corresponding  set  of 
uniform  spatial  samples  which  have  spectral  coeffi- 
cients that  approximate  the  uniform  spectrum  of  the 
delta  function. 

The  procedure  in  one  dimension  is  as  follows;  Let 

m-1 

6(x-Xo)=  &(,sin^x  with  0-^x^tt',  (26) 

f=i 


(29) 


More  explicity,  if  we  arbitrarily  select  the  number 
,of  nonzero  coefficients  in  the  sine  series  (w-  1)  to  be 
even,  an  expression  for  the  sample  coefficients  as  a 
function  of  the  fraction  of  the  center  interval  at  which 
Xg  is  located  is 


which  is  readily  solved  for  the  coefficients 

bp=(2/T!)sin(pXg)  . (27) 

An  inverse  discrete  Fourier  transform  then  gives  a set 
of  sample  coefficients 

E sin/)XoSin^/)a;^^  , a = 1,  2,  . . . , m - 1 . 

(28) 

This  set  of  coefficients  wili  not  be  exactly  normalized 
for  all  values  of  Xq  because  of  truncation  errors,  so  a 
better  approximation  to  the  delta  function  is  achieved 
by  the  set  of,  coefficients 


! m — \ \1 

p — 

— H — + € 1 

m 

V 2 j\ 

(30) 


This  direct  expansion  of  the  delta  function  exhibits 
the  troublesome  Gibbs  oscillations  in  its  power  spec- 
trum as  a consequence  of  the  truncation  of  the  infinite 
series.  Lanczos®  recommends  what  he  terms  “a 
smoothing”  which  applies  a sin(x)/x  weighting  function 
to  the  coefficients  of  the  truncated  trigonometric  series. 
These  smoothed  coefficients  are  given  by 


g'ji)  = gji) 


sin[T7(ffl  - zm)/zm] 
u(o!  - 


(31) 


TABLE  IV.  Unsmoothed  coefficients  for  an  8-point  delta-function  expansion. 


e 

1 

2- 

3 

a 

4 5 

6 

7 

8 

C.  G, 

0 

0.0000 

0,0000 

0.0000 

1,0000 

0.0000 

0,0000 

0.0000 

0.0000 

4.0000 

0.1 

-0.0148 

0,0354 

-0.0829 

0,9867 

0. 1087 

-0,0468 

0.0245 

-0.0109 

4.1013 

0.2 

-0. 0271 

0,0642 

-0.1447 

0.9412 

0,2339 

-0.0946 

0.0488 

-0.0216 

4. 2018 

0.3 

- 0.  0359 

0.0843 

-0.1838 

0.8655 

0.3694 

-0.1388 

0.0702 

-0.0309 

4.3017 

0.4 

-0.0407 

0.0946 

-0.2004 

0.7639 

0.5082 

- 0. 1743 

0.0864 

- 0.  0378 

4.  4010 

0.5 

- 0.  0412 

0.0950 

-0, 1962 

0,6424 

0. 6424 

- 0. 1962 

0.0950 

-0.0412 

4. 5000 

0.6 

-0.0378 

0.0864 

- 0. 1742 

0.5082 

0.7639 

-0. 2004 

0.0946 

-^0.0407 

4.5990 

0.7 

- 0.  0309 

0.0702 

-0.1388 

0.3694 

0. 8655 

-0.1838 

0 . 0843 

-0.0359 

4.  6983 

0.8 

-0.0216 

0.0488 

-0.0946 

0.2339 

0.9412 

-0.1447 

0.  0642 

-0.0271 

4.7982 

0.9 

-0.0109 

0.0245 

-0.0468 

0,1087 

0.9867 

-0.0829 

0.0354 

-0.0148 

4. 8987 

1.0 

0,0000 

0.0000 

0.0000 

0.0000 

1. 0000 

0.0000 

0.0000 

0. 0000 

5.0000 
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TABLE  V.  Smoothed  coefficients  for  an  8-point  delta-function  expansion. 


e 

1 

2 

3 

4 

a 

5 

6 

7 

8 

C.  G, 

0 

0.0000 

0.0000 

0.0000 

1.0000 

0.0000 

0.0000 

0.0000 

0.0000 

4.  0000 

0. 1 

-0.0039 

0.0201 

-0.0690 

0. 9735 

0.1073 

-0.0390 

0.0139 

-0.0029 

4. 1001 

0.2 

-0.0071 

0.0361 

-0.1191 

0.  9182 

0.2281 

-0.0779 

0.0274 

-0.0057 

4.2001 

0.3 

-0.0093 

0.0470 

-0.1501 

0.  8373 

0.3574 

- 0. 1133 

0.0391 

-0.0080 

4.3001 

0.4 

-0,0105 

0.0524 

-0.1628 

0.7352 

0.4891 

- 0. 1415 

0.0479 

-0.0098 

4.4000 

0.5 

-0.0106 

0.0526 

-0. 1591 

0.6171 

0.6171 

-0.1591 

0.0526 

-0.0106 

4.5000 

0.6 

-0.0098 

0.0479 

- 0. 1415 

0.4891 

0.7352 

-0.1628 

0.0524 

-0.0105 

4.6000 

0.7 

-0.0080 

0.0391 

-0.1133 

0.3574 

0.8373 

-0.1501 

0.0470 

-0.0093 

4.6999 

0.8 

-0.0057 

0.0274 

-0,0779 

0.2281 

0.9182 

-0.1191 

0.0361 

-0,0071 

4.  7999 

0.9 

-0.0029 

0.0139 

-0.0390 

0.1073 

0.9735 

-0.0690 

0.0201 

-0,0039 

4.8999 

1.0  . 

0. JOOO 

0.  0000 

0.0000 

0.0000 

1.0000 

0.0000 

0.0000 

0.0000 

5.0000 

with  the  corresponding  normalized  coefficients, 

^'«(€)  = ^„(e)/ E . (32) 

' a=l  , 

Sets  of  unsmoothed  and  smoothed  coefficients  for  an 
8-sample  expansion  are  tabulated  in  Tables  IV  and  V. 

An  indication  of  the  accuracy  with  which  these  expan- 
sions represent  the  delta  function  is  given  by  the  center 
of  gravity  of  the  set,  also  listed  in  Tables  IV  and  V, 
which  is  a measure  of  the  phase  error  in  the  low-fre- 
quency region  of  the  transformed  spectrum . The  ac- 
curacy with  which  the  center  of  gravity  tracks  the  actual 
location  of  the  delta  function  is  considerably  better  in 
the  smoothed  sets  (±0.0001  of  an  interval),  than  in  the 
unsmoothed  sets  (±0.0018  of  an  interval).  However, 
the  maximum  error  for  either  has  a negligible  effect 
on  the  computed  pattern  as  given  by  Eq.  (8). 

A look  at  the  spectrum  of  the  sample  coefficients 
gives  an  even  better  understanding  of  the  fidelity  of 
representation.  The  maximum  deviation  of  the  power' 
spectrum  from  the  uniform  spectrum  of  the  delta  func- 
tion occurs  for  the  midinterval  set.  Relative  power 
spectra  for  several  unsmoothed  midinterval  sam- 
ple sets  are  plotted  in  Fig.  7;  corresponding  spec- 
tra for  smoothed  sets  are  shown  in  Fig.  8.  The  ex- 
pansion, in  effect,  replaces  the  element  with  a direc- 
tional subarray  whose  directivity  pattern  is  represented 
by  the  power  spectrum  of  the  sample  set.  This  varies 
from  a uniform  spectrum  when  the  location  of  the  ele- 
ment coincides  with  a sample  position  to  the  maximum 
deviation  shown  in  Figs.  7 and  8 when  the  location  of  the 
element  falls  in  the  middle  of  a sampling  interval. 

Qualitatively,  the  errors  indicated  by  the  extreme 
deviations  of  the  power  spectrum  can  be  related  to  the 
array  pattern  error  of  Eq.  (8).  For  example,  the  un- 
smoothed 4-point  expansion  shows  a + 1.2-dB  excursion 
at  the  midpoint  of  the  transform  spectrum.  If  this 
change  in  amplitude  from  1.0  to  1. 14  were  considered 
to  be  distributed  uniformly  in  a random  manner  over 
the  element  set,  it  would  give  rise  to  a CT„  = -n  (0. 14)^ 

= 0.0016.  Thus,  even  this  crude  4-point  expansion 
produces  an  error  nearly  30  dB  below  the  mean  side- 
lobe  level.  Because  the  representation  is  so  good,  the 
dominant  critierion  for  selection  of  the  size  of  the  set 
for  expansion  of  the  delta  function  is  that  of  the  usable 
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fraction  of  the  transform  field.  Figure  9 summarizes 
the  pass  bands  of  the  spectra  of  Figs.  7 and  8.  The 
unsmoothed  expansions  provide  a significantly  larger 
usable  band,  and  the  smoothed  versions  would  only  be 
required  for  very  precise  computations  or  for  super 
gain  arrays  with  large  positive  and  negative  weights, 
where  the  array  sum  becomes  very  sensitive  to  errors 
in  either  phase  or  amplitude. 

Extending  the  expansion  to  two  or  three  dimensions  is 
straightforward.  For  two  dimensions 

771-1  771-1 

5{x-Xj,  y-y^)=E  E ^o®sin(/)A;)sin(<7y)  , 

#»=1  G=1 


FIG.  7.  Transform  spectra  for  unsmoothed  midinterval  delta- 
function  expansion. 
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Fraction  of  transform  field 

FIG.  8.  Transform  spectra  for  t smoothed  midinterval  delta- 
function  expansion. 


with  0 77,  0<y>TT, 


m 


Wt**  1.  1 j ^ \ ™ . V 


7 ffi-i 


FIG.  9,  Usable  fraction  of  the  transform  field  versus  the  num- 
ber of  coefficients  in  the  delta-function  expansion. 


the  first  step  in  applying  the  DFT  to  the  problem  is  to 
select  the  number  of  samples  for  the  transform  field. 
Figure  10  illustrates  the  relationship  between  the  trans- 
form parameters.  In  general,  the  maximum  range  of 
^max>  '''ill  1^®  larger  than  the  visible  range,  4/\,  by 
a factor  i/s  1 and  the  maximum  range  in  the  spatial 
coordinate,  will  be  larger  than  the  aperture  di- 

mension X>  by  a factor  ya  1.  The  sampling  interval  in 
X is  ^ = ^'''1 11'®  sampling  interval  in  k is  C^k 

=^max-  Combining  these  relationships,  the  minimum 
number  of  samples  in  the  DFT  field  is  given  by  L„,i„ 

= AyvD/\. 

The  selection  of  y is  predicted  on  two  criteria: 

(1)  y is  greater  than  or  equal  to  the  sampling  fre- 
quency ratio  selected  for  interpolation  of 

(2)  y 3;  (D  + m£^)/D,  where  m is  the  number  of  coef- 


If  all  of  the  elements  of  the  array  have  uniform  ampli- 
tude weights,  there  may  be  a computational  advantage  to 
precomputing  and  storing  this  four-dimensional  matrix 
in  a look-up  table,  providing  the  required  memory  is 
available.  Otherwise,  the  sample  coefficients  can  be 
readily  generated  from  a much  smaller  two-dimen- 
sional table  corresponding  to  Table  IV  or  V. 

Computing  time  for  this  expansion  is  proportional 
to  N{m-  if  for  a two-dimensional  array  so  the  selec- 
tion of  m depends  on  the  trade  off  between  the  increase 
in  the  computing  time  for  generating  the  expansion  co- 
efficients as  m is  increased  and  the  reduction  in  DFT 
computing  time  achieved  by  the  increase  in  the  fraction 
of  the  transformed  field  which  is  usable. 


IX.  APPLYING  THE  DIGITAL  FOURIER  TRANSFORM 

Now  that  criteria  have  been  established  for  uniformly 
sampling  both  the  k domain  and  the  array  aperture. 


A 


DIGITAL 

FOURIER 

TRANSFORM 


MINIMUM  NUMBER  OF  POINTS  IN  FFT 
, " MAX  _ " MAX  4 V y 0 


A K 


Ax 


ARRAY  APERTURE 


V J 

^ ''  MAX 

1*4 ^ — ARRAY  APERTURE  * 

1 M'N  1 1 1 1 1 1 1 1 

= 0 

1 1 1 1 1 1 1 

^Ax=  ' = ^ 

FIG.  10.  Parameter  relationships  in  the  discrete  Fourier 
transform. 
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FIG.  11.  Computation  times  for  various  algorithms  versus  the 
number  of  samples  in  the  output  field. 


ficients  used  for  the  5{xj)  expansion. 

Of  these  two  criteria,  usually  (1)  is  the  controlling  one. 

A single  criterion  is  used  for  the  selection  of  p; 

4/X  + wAfe 

> 

usable  fraction  of  for  selected  m ’ 

where  n is  the  number  of  samples  used  for  inter- 
polation. These  above  criteria  apply  separately  to  the 
separable  transformations.  That  is,  A„min, 

are  related  to  each  other  but  independent  of  i^min, 
p,,  y^,  DyOrL^min,  v„  y^,  D^. 

In  order  to  exploit  the  use  of  packaged  FFT  algorithms 
of  radix  2,  it  is  necessary  to  select  an  integer  power  of 
2 field  size;  i.e.,  A = 2'’>Lmin.  The  scale  factor  Ajc 
is  used  to  relate  array-element  locations  to  sample  in- 
dices and  a fractional  interval  for  delta  function  ex- 
pansion. Note  that  the  scale  factors  for  x,  y,  z will 
usually  be  different  except  for  the  special  case  when 
the  aperture  size  is  the  same  in  two  or  more  coordinates. 

The  advantages  of  using  a symmetrical  array  distri- 
bution appear  in  the  application  of  the  FFT  algorithm 
also.  By  imposing  the  condition  of  a real  symmetrical 
distribution  on  the  array  aperture,  a cosine  transform 
which  utilizes  a complex  FFT  working  on  i as  many 
sample  points  with  a corresponding  savings  in  time  of 
about  4 ; 1 may  be  used.  An  algorithm  for  the  cosine 
transform  is  described  by  Cooley  et  al.  ’ This  savings 
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of  4 ; 1 becomes  16;  1 for  a two-dimensional  array  and 
64;  1 for  a three-dimensional  array  and  thus  is  well 
worth  consideration.  If  the  cosine  transform  is  used, 
only  one  quadrant  of  a two-dimensional  array  aperture 
would  be  used  and  thus  in  the  expansion  of  elements  in- 
to the  uniform  sample  set,  coefficients  of  elements  in 
the  (+,  +)  quadrant  would  overlap  into  adjacent  quad- 
rants. These  coefficient  values  are  merely  folded  back 
into  the  (+,  +)  quadrant  to  complete  the  array-distri- 
bution matrix. 

X.  COMPUTATION-TIME  GUIDELINES 

A simple  computer  model  having  a memory  cycle 
time  of  2.5  psec,  addition/subtraction  of  6 psec, 
multiplication/division  of  7 psec,  and  a cosine  routine 
of  135  psec  will  be  assumed  for  the  purpose  of  making 
a relative  comparison  of  computation  times  for  the 
techniques  which  have  been  discussed  above. 

Figure  11  displays  approximate  computation  times 
for  a two-dimensional  array  for  the  interpolation 
algorithm,  FFT,  and  direct  computations  as  a function 
of  logioM,  the  number  of  samples  in  the  field. 

One  set  of  lines  is  for  the  direct  computation  with  N, 
the  number  of  elements,  as  a parameter  assuming 
/=  16,  a 16-point  interpolation.  The  use  of  a stored  • 
cosine  table  along  with  array  symmetry  is  assumed  for 
these  times.  ‘ The  other  set  of  lines  represent  inter- 
polation times  for  2-  16  points  of  interpolation  filtering. 
The  complex  and  cosine  FFT  times  are,  represented  as 
bands  signifying  that  the  FFT  size  changes  in  discrete 
steps  of  a power  of  two  in  both  dimensions  and  times 
are  bounded  by  the  bands  shown.  The  FFT  times/do 
not  include  the  increased  computation  required  by  the 
factor  in  comparison  with  the  direct  computation. 

Not  displayed  in  Fig.  11  are  the  times  for  the  delta- 
function  expansions  required  for  the  FFT  computations. 
For  the  two-dimensional  case  the  time  will  be  given 
approximately  by  N(m-  l)^x25  psec  where  (m-  1)  is 
the  number  of  coefficients  in  the  expansion.  For  filled 
arrays  (average  spacing  ix)  this  time  is  comparable  to 
the  cosine  FFT  time  at  the  Nyquist  sampling  rate  for 
w - 1 = 8.  As  such  it  represents  a minor  but  not  neg- 
ligible part  of  the  overall  computation  time. 

The  sequential  use  of  the  interpolation  algorithms 
will  be  illustrated  in  a specific  example. 

XI.  COMPUTATION-TIME  EXAMPLE 

To  summarize  the  reduction  in  computing  time  which 
can  be  effected  by  the  techniques  that  have  been  dis- 
cussed, let  us  choose  a hypothetical  array  for  which 
computing  time  would  be  significant. 

Take  as  an  example  a set  of  N=400  uniformly 
weighted  elements  randomly  distributed  over  a plane 
aperture  of  20XX20X.  Our  goal  will  be  to  represent 
the  array  factor  by  a data  array  covering  the  two- 
dimensional  visible  region  (a  fully  steered  hemisphere) 
to  an  accuracy  of  3%  (-  30  dB).  As  a basis  for  esti- 
mating computing  time  we  will  use  the  cycle  times  as- 
sumed in  the  previous  section. 
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The  total  visible  range  of  k is  4/\.  Therefore,  the 
total  number  of  Nyquist  intervals  in  one  dimension  of 
our  20X  array  is  80. 

From  Table  II  we  find  that,  for  a data  array  (single- 
point interpolation)  the  sampling  frequency  must  be 
a 13.  SxNyquist  to  obtain  an  error  <-30  dB.  Thus,  the 
number  of  elements  in  one  row  of  our  data  matrix  is 
80X  13.8=  1104  for  a total  set  of  M=  (1104)^=  1.  22x  10®. 

A.  Direct 

The  direct  computation  of  Eq.  (3)  then  requires  2MN 
X (computation  time  per  operation)  = 2x 400x  1.  22x  10® 
x(2xl35+7x2,5+2x7  + 4x6)x  10‘®  = 1.  58  X 10®  sec  = 44  h. 
One  might  hesitate  before  embarking  on  such  a com- 
puting exercise! 

B.  Symmetry 

The  first  step  is  to  introduce  symmetry  into  the  array 
distribution.  The  time  now  becomes  (M/4)x(iVf/4)x(2 
xl35+5x2.5  + 4x7+6)  = i(1.22xl0®)x|(400)x(316) 

= 9.  64x  10®  sec  = 2.  7 h.  This  is  a good  start  at  improv- 
ing the  computer  time. 

C.  Look-up  table 

Now  use  a stored  table  instead  of  the  cosine  routine. 
For  an  error  of  - 30  dB  a quantization  step  of  0.  1 rad 
is  required.  The  range  of  the  argument  of  the  function 
for  the  symmetrical  case  is  2itD/\.  Thus,  a table  of 
1250  elements  will  be  required  at  an  initial  set-up  time 
of  1250X  135 -t-  2.  5-t-  6x  10‘®  = 0. 18  sec— negligible. 

The  computation  time  may  be  obtained  directly  from 
Fig.  11.  N=400,  /=  1.  22  X 10®  gives  a time  of  10®*® 

= 1.  12x10®  psec  = 0.31  h. 

D.  Direct  with  interpolation 

Finally,  introduce  interpolation  using  a 16,  4,  2,  2 
sequence.  Sampling  the  final  field  at  Nyquist  rate  will 
require  iff  = 80  x 80  = 6400  samples,  1. 06  Nyquist  would 
be  6784.  From  Fig.  11,  the  time  required  for  400 
elements  will  be  fi  = 9.  55  sec  indicated  by  point  1 on  the 
figure.  The  first  interpolation  step,  a 16-point  weighting 
window,  at  point  2,  required  lo®*®®  p,sec,  f2  = 3.63  sec. 
The  second  step,  a 4-point  window  at  point  3 takes  fg 
= 3.  02  sec.  Next,  a 2-point  window  at  4 requires 
= 3.  39  sec.  Finally,  for  the  last  step,  also  a 2-point 
window,  at  5 12.9  sec.  The  total  time  fi-hfg  + fg+fi 

-1-^5  = 32.  5 sec. 

E.  Cosine  FFT  with  interpolation 

An  alternate  route  is  fhe  FFT  approach.  We  select 
an  8-coefficient  unsmoothed  expansion  for  the  element 
locations  for  which,  from  Fig.  9,  a of  1.  25  is  ap- 
propriate. 

After  some  trial  and  error  with  the  minimum  y ratios 
of  Table  II  and  the  relative  computing  times  involved 
in  the  FFT  and  interpolation  steps,  a sequence  of 
cosine  FFT,  and  a 4-point  interpolation  followed  by  two 


2-point  interpolations  is  selected  with  the  controlling 
7=3.63  being  the  required  sampling  rate  at  the  output 
of  the  4-point  interpolation.  The  sequence  is  shown 
as  steps  I,  II,  III,  and  IV  in  Fig.  ,11.  This  selection 
gives  precisely  the  required  y at  steps  IT  and  III  and 
provides  the  minimum  value  achievable  at  step  IV, 
slightly  in  excess  of  that  required  but  less  than  that  re- 
sulting from  the  above  direct  method.  The  value  of  y 
for  the  cosine  FFT  is  3.63/2  = 1.  81.  From  these  values 
of  V and  y,  I,min  = 4x  l.  25x  1.  81  x 20  = 181  therefore, 

L = 253.  The  cosine  transform  is  carried  out  with  i 
as  many  samples,  so  the  cosine  transform  involves  a 
64-point  FFT  for  one  dimension.  For  the  two-dimen- 
sional transform  the  total  number  of  points  is  64® 

= 4096.  Processing  time  will  be  4096xlog2(4098)x  84 
psec  = 4.  13  sec.  The  time  required  for  the  8-coefficient 
expansion  for  400  elements  is  0.  64  sec.  The  4-point 
interpolation  of  step  II  takes  2.2  sec,  step  III,  2.5  sec, 
and  step  IV,  10  sec,  for  a total  of  19.  5 sec. 

This  represents  a savings  of  about  40%  over  the  di- 
rect/interpolation method  but  at  the  cost  of  somewhat 
greater  programming  complexity.  For  arrays  with 
more  elements  the  savings  would,  of  course,  be 
greater.  If  the  complex  form  of  the  FFT  were  used 
directly  instead  of  the  more  efficient  cosine  FFT 
algorithm  the  computing  time  would  be  essentially  the 
same  as  for  the  direct/interpolation  method  in  this  ex- 
ample. 

It  is  interesting  to  note  from  Fig.  11  that  the  use  of 
interpolation  offers  more  than  an  order-of-magnitude 
decrease  in  computing  time  for  the  y required  for  ar- 
ray pattern  plotting  or  for  numerical  representation 
over  even  the  very  efficient  cosine  FFT  process. 

It  is  also  apparent  that  either  the  32.  5-  or  19.5-sec 
times  are  much  more  tolerable  figures  than  any  of' the 
other  computing  times,  particularly  more  so  than  the 
44-h  “brute-force”  approach,  a ratio  of  8000;  1. 

XII.  CONCLUSION 

In  the  face  of  a savings  of  8000 : 1 what  else  can  one 
conclude  but  that  efficient  transform  algorithms  com- 
bined with  interpolation  algorithms  really  pay  their  way 
in  array-pattern  computation. 


*R.  E,  Collin  and  X.  Zucker,  in  Antenna  Theory  Part  edit- 
ed by  R.  E.  Collin  (McGraw-Hill,  New  York,  1969),  p.  141. 

®R.  E.  Collin,  Ref.  1.  p.  147. 

®R.  E.  Collin,  Ref.  1.  231,  Eq.  (6.40). 

^R.  B.  Blackman  and  J.  W.  Tukey,  “The  Measurement  of 
Power  Spectra,”  Bell  Syst."' Tech.  J.  37,  76  (1958). 

®C.  Lanczos,  Applied  Analysis  (Prentice-Hall,  Englewood 
Cliffs,  1956),  pp.  229-241. 

®C.  Lanczos,  Ref.  5,  p,  225. 

^J.  W.  Cooley,  P.  A.  Lewis,  and  P.  D.  Welch,  “The  Fast 
Fourier  Transform  Algorithm:  Programming  Considera- 
tions in  the  Calculation  of  Sine,  Cosine  and  Laplace  Trans- 
forms,” in  Digital  Signal  Processing,  edited  by  L.  R.  Rabi- 
ner  and  C.  M.  Rader  (IEEE,  New  York,  1972). 


J.  Acoust.  Soc.  Am.,  Vol.  61,  No.  3,  March  1977 


DISTRIBUTION  LIST 


Chief  of  Naval  Research 
Department  of  the  Navy 
Arlington,  Virginia  22217 


Code 

200 

(1) 

Code 

222 

(2) 

Code 

102IP 

(1) 

Code 

1020SC 

(1) 

Code 

460 

(1) 

Code 

480 

(1) 

Code 

481 

(1) 

Code 

486 

(1) 

Director 

Office  of  Naval  Research 
Branch  Office 
1030  East  Green  Street 

Pasadena,  Calfornia  91101  (1) 


Commander 

Naval  Sea  Systems  Command 


Washington,  D.C. 

20362 

Code  03E 

(1) 

Code  034 

(1) 

Code  0342 

(1) 

Code  036 

(1) 

Code  06H1 

(1) 

Code  06H2 

(1) 

Code  09G3 

(1) 

Code  92 

(1) 

Code  662C14 

(1) 

PMS  395-4 

(1) 

PMS  402  B 

(1) 

Code  03423  (Mr. 

Francis  J.  Romano) 

(1) 

Commander 

Naval  Air  Systems 

Command 

Washington,  D.C. 

20361 

Code  370 

(1) 

Code  264 

(1) 

Commander 
Naval  Electronics 
Systems  Command 
Washington,  D.C.  20360 

Code  PME-124  (1) 

Chief  of  Naval  Material 
Department  of  the  Navy 
Washington,  D.C.  20360 

Code  PM-4  (1) 

Code  034  (1) 

Code  ASW121  (1) 

Chief  of  Naval  Operations 
Department  of  the  Navy 
Washington,  D.C.  20350 

Code  Op  03  (1) 

Code  Op  32  (1) 

Code  Op  098  (1) 

Code  Op  02  (1) 

Code  Op  095  (1) 

Code  Op  23  (1) 

Code  Op  967  (1) 


Commanding  Officer 
Naval  Ocean  Research  8 

Development  Activity  (NORDA) 

National  Space  Technology  Laboratories 
Bay  St.  Louis,  Mississippi  39529 


Code  100  (1) 

Code  460  (1) 

Code  500  (1) 

Oceanographer  of  the  Navy 
The  Hoffman  Building 
200  Stoval  Street 

Alexandria,  Virginia  22332  (1) 


Director  Strategic  Systems 
Projects  Office  (PM-1) 

Department  of  the  Navy 
Washington,  D.C.  20360 

Code  NSP-20  (1) 

Director 

Defense  Research  5 Engineering 
The  Pentagon 
Washington,  D.C.  20301 
Assistant  Director  (Sea 


Warfare  Systems)  (1) 

Assistant  Secretary  of  the  Navy 
(Research  § Development) 

Department  of  the  Navy 

Washington,  D.C.  20350  (1) 

U.S.  Naval  Oceanographic  Office 
Washington,  D.C.  20373 

Code  1640  (1) 

Code  3440  (1) 

Commander  Operational  Test  8 
Evaluation  Force 
U.S.  Naval  Base 

Norfolk,  Virginia  23511  (1) 

Commander,  Submarine  Force 
U.S.  Pacific  Fleet 
Fleet  Post  Office 

San  Francisco,  California  96601  (1) 

Commander 

Submarine  Group  FIVE 

Fleet  Station  Post  Office 

San  Diego,  California  92132  (1) 

Commander 
Third  Fleet 

U.S.  Pacific  Fleet,  FPO 

San  Francisco,  California  96610  (1) 

Commander 

Submarine  Development  Group  ONE 
Fleet  Post  Office 

San  Diego,  California  92132  (1) 

Commander 

Submarine  Development  Group  TWO 
Naval  Submarine  Base-New  London 
Groton,  Connecticut  06340  (1) 


I 


DISTRIBUTION  LIST  (Continued) 


Commander,  Surface  Force 
U.S.  Atlantic  Fleet 

Norfolk,  Virginia  235II  (1) 

Commander,  Surface  Force 
U.S.  Pacific  Fleet 

San  Diego,  California  92155  (I) 

Reprint  Custodian 

Department  of  Nautical  Science 

U.S.  Merchant  Marine  Academy 

Kings  Point,  New  York  11024  (1) 

Deputy  Commander 
Operational  Test  § 

Evaluation  Force,  Pacific 
U.S.  Naval  Air  Station 


San  Diego,  California  92135  (1) 

Commander 

Naval  Ship  Reserach  5 
Development  Center 

Bethesda,  Maryland  20084  (1) 

Naval  Civil  Engineering  Laboratory 
Port  Hueneme,  California  93041 

Code  L40  (1) 

Code  L42  (1) 

Naval  Facilities 

Engineering  Command 
Washington,  D.C,  20390 

Code  03  (1) 

Code  032C  (1) 

Commanding  Officer 
U.S.  Naval  Air 

Development  Center 

Warminster,  Pennsylvania  18974  (1) 

Commander 

Naval  Ocean  Systems  Center 
San  Diego,  California  92152 

Code  6700  (2) 

Officer  In  Charge 
Naval  Ship  Research  § 

Development  Center 

Annapolis,  Maryland  21402  (1) 

Commanding  Officer,  Naval 
Coastal  Systems  Laboratory 
Panama  City,  Florida  32401  (1) 


Commander 
Naval  Surface 

Combat  Systems  Center 
White  Oak 

Silver  Spring,  Maryland  20910  (1) 


Director 

U.S.  Naval  Research  Laboratory 
Washington,  D.C.  20375 


Code 

2620 

CD 

Code 

2627 

CD 

Code 

8000 

CD 

Code 

8100 

CD 

Commanding  Officer 

Naval  Underwater  Systems  Center 

Newpor,t,  Rhode,  Island  20844 
Attn:  John  D ' A I bo  ra 


Commanding  Officer 

Naval  Underwater  Systems  Center 

New  London,  Connecticut  06320 

Code  900  (1) 

Code  905  (1) 

Code  910  (1) 

Code  930  (1) 

Code  960  (1) 

Commanding  Officer 

Naval  Training  Equipment  Center 

Orlando,  Florida  32813 

Tech  Library  (1) 


Chief  Scientist 
Navy  Underwater 

Sound  Reference  Division 
U.S.  Naval  Research  Laboratory 
P.O.  Box  8337 
Orlando,  Florida  32806 

Superintendent 

U.S.  Naval  Postgraduate  School 
Monterey,  California  93940  (1) 

Director 

Defense  Documentation  Center 
(TIMA) , Cameron  Station 
5010  Duke  Street 

Alexandria,  Virginia  22314  (12) 


Executive  Secretary 
National  Academy  of  Sciences 
2101  Constitution  Avenue,  N.W. 
Washington,  D.C.  20418  (1) 

Supreme  Allied  Commander 
U.S.  Atlantic  Fleet 
ASW  Research  Center 
APO 

New  York,  New  York  09019 
Via:  ONR  210 

CNO  0P092D1 

Secretariat  of  Military 

Information  Control 

Committee  (1) 


II 


DISTRIBUTION  LIST  (Continued) 


National  Oceanic  & Atmospheric 
Administration 
6001  Executive  Boulevard 
Rockville,  Maryland  20852  (1) 

Director  of  Naval  Warfare  Analysis 
Institute  of  Naval  Studies 
1401  Wilson  Boulevard 

Arlington,  Virginia  22209  (1) 

Institute  for  Defense  Analyses 
400  Army-Navy  Drive 

Arlington,  Virginia  22202  (1) 

Director 
Woods  Hole 

Oceanographic  Institution 
Woods  Hole,  Massachusetts  02543  (1) 

Meterological  6 

Astrophysical  Abstracts 
310  E.  Capitol  Street 

Washington,  D.C.  20003  (1) 

Director 

Applied  Physics  Laboratory 
University  of  Washington 
1013  East  40th  Street 

Seattle,  Washington  98105  (1) 

National  Science  Foundation 
Washington,  D.C.  20550  (1) 

Director 
Lamont-Doherty 

Geological  Observatory 
Torrey  Cliff 

Palisades,  New  York  10964 


Director 

Applied  Research  Laboratory 
Pennsylvania  State  University 
P.O.  Box  30 

State  College,  Pennsylvania  16802  (1) 


Director 

University  of  Texas 
Applied  Research  Laboratory 
P.O.  Box  8029 

Austin,  Texas  78712  (1) 


TACTEC 

Battel le  Columbus  Laboratories 
505  King  Avenue 

Columbus,  Ohio  43201  (1) 


Director 
Institute  of 

Ocean  Science  8 Engineering 
Catholic  University  of  America 
Washington,  D.C.  20017  (1) 


Director 

Marine  Research  Laboratories 

c/o  Marine  Studies  Center 

University  of  Wisconsin 

Madison,  Wisconsin  53706  (1) 


Office  of  Naval  Research 
Resident  Representative 
c/o  University  of  California, 
San  Diego 


La  Jolla,  California  92093  (1) 


University  of  California, 

San  Diego 

La  Jolla,  California  92093 

MPL  Branch  Office  (5) 


Director 

College  of  Engineering 
Department  of  Ocean  Engineering 
Florida  Atlantic  University 
Boca  Raton,  Florida  33431 


ONR/MPL  General  May  1977 


III 


Marine  Physical  Laboratory 
MPL-U-23/76 

EFFICIENT  COMPUTATION  OF  ARRAY  PATTERNS  by  Victor  C. 
Anderson,  University  of  California,  San  Diego,  Marine 
Physical  Laboratory  of  the  Scripps  Institution  of 
Oceanography,  San  Diego,  California  92152, 
Journal  of  the  Acoustical  Society  of  America,  Vol. 

61,  No.  3,  March  1977, 

The  impact  of  a symmetrical  array  geometry,  the  use  of  a 
quantized  stored  cosine  function,  the  exploitation  of 
digital  Fourier  transform  algorithms,  and  the  application 
of  trigonometric  interpolation  in  the  computation  of  array 
patterns  is  discussed.  Careful  selection  of  parameters 
permits  sampling  the  array  pattern  only  6%  above  the 
theoretical  Nyquist  limit.  Reconstruction  of  array 
patterns  showing  -20,  -30,  and  -40-dB  relative 
interpolation  errors  are  presented.  A saving  of  8000:1  in 
computation  time  over  direct  "brute-force"  array-pattern 
computation  is  illustrated  for  a hypothetical  array. 


II.  Signal  Processing 
a.  Instrumentation 


Victor  C,  Anderson 


1 Sponsored  by 
Office  of  Naval  Research 
N00014-7S-C-0749 
and 

Advanced  Research  Projects 
Agency 

N00014-76-C-1080 


Marine  Physical  Laboratory 
MPL-U-23/76 

EFFICIENT  COMPUTATION  OF  ARRAY  PATTERNS  by  Victor  C. 
Anderson,  University  of  California,  San  Diego,  Marine 
Physical  Laboratory  of  the  Scripps  Institution  of 
Oceanography,  San  Diego,  California  92152, 
Journal  of  the  Acoustical  Society  of  America,  Vol. 

61,  No.  3,  March  1977. 

The  impact  of  a symmetrical  array  geometry,  the  use  of  a 
quantized  stored  cosine  function,  the  exploitation  of 
digital  Fourier  transform  algorithms,  and  the  application 
of  trigonometric  interpolation  in  the  computation  of  array 
patterns  is  discussed.  Careful  selectiotJ  of  parameters 
permits  sampling  the  array  pattern  only  6%  above  the 
theoretical  Nyquist  limit.  Reconstruction  of  array 
patterns  showing  -20,  -30,  and  -40-dB  relative 
interpolation  errors  are  presented.  A saving  of  8000; 1 in 
computation  time  over  direct  "brute-force"  array-pattern 
conq)utation  is  illustrated  for  a hypothetical  array. 


UNCLASSIFIED 


II.  Signal  Processing 
a.  Instr‘Amentation 


Victor  C.  Anderson 


Sponsored  by 
Office  of  Naval  Research 
N00014-75-C-0749 
and 

Advanced  Research  Projects 
Agency 

N00014-76-C-1080 


UNCLASSIFIED 


Marine  Physical  L.nboratory 
MPL-U-23/76 

II.  Signal  Processing 
a.  Instrumentation 

EFFICIt.'T  COMPUTATION  OF  ARRAY  PATTIiRNS  by  Victor  C. 
Anderson,  University  of  California,  San  Diego,  Marine 
Physical  Laboratory  of  the  Scripps  Institution  of 
Oceanography,  San  Diego,  California  92152. 
Journal  of  the  Acoustical  Society  of  America,  Vol. 
61,  No,  3,  March  1977, 

Victor  C.  Anderson 

The  impact  of  a symmetrical  array  geometry,  the  use  of  a 
quantized  stored  cosine  function,  t!.e  exploitation  of 
digital  Fourier’  transform  algorithms,  and  the  application 
of  trigonometric  interpolation  in  the  computation  of  array 
patterns  is  discussed.  Careful  selection  of  parameters 
permits  sampling  the  array  pattern  only  6%  above  the 
theoretical  Nyquist  limit.  Reconstruction  of  array 
patterns  showing  -20,  -30,  and  -40-dB  relative 
interpolation  errors  are  presented,  A saving  of  8000:1  in 
computation  time  over  direct  "brute-force"  array-pattern 
computation  is  illustrated  for  a hypothetical  array. 

Sponsored  by 
Office  of  Naval  Research 
N00014-7S-C-0749 
and 

Advanced  Research  Projects 
Agency 

N00014-76-C-1080 

UNCLASSIFIED 

Marine  Physical  Laboratory 
MPL-U-23/76 

II.  Signal  Processing 
a.  Instrumentation 

EFFICIENT  COMPUTATION  OF  ARRAY  PATTERNS  by  Victor  C. 
Anderson,  University  of  California,  San  Diego,  Marine 
Physical  Laboratory  of  the  Scripps  Institution  of 
Oceanography,  San  Diego,  California  92,152, 
Journal  of  the  Acoustical  Society  of  America,  Vol. 
Cl,  No.  3,  March  1977. 

Victor  C.  Anderson 

The  impact  of  a symmetrical  array  geometry,  the  use  of  a 
quantized  stored  cosine  function,  the  exploitation  of 
digital  Fouiier  transform  algorithms,  and  the  application 
of  trigonometric  interpolation  in  the  computation  of  array 
patterns  is  discussed.  Careful  selection  of  parameters 
permits  sampling  the  array  pattern  only  6%  above  the 
theoretical  Nyquist  limit.  Reconstruction  of  array 
patterns  showing  -20,  -30,  and  -40-dB  relative 
interpolation  errors  are  presented.  A saving  of  8000:1  in 
computation  time  over  direct  "brute-force"  array-pattern 
computation  is  illustrated  for  a hypothetical  array. 

Sponsored  by 
Office  of  Naval  Research 
N00014-75-C-0749 
and 

Advanced  Research  Projects 
Agency 

N00014-76-C-1080 

UNCUSS  IFIED 


